UNITS {
  (mV) = (millivolt)
  (mA) = (milliamp)
  (S) = (siemens)
}

NEURON {
  SUFFIX mnsoma43
  USEION Ca READ eCa WRITE iCa VALENCE 1
  USEION Na READ eNa WRITE iNa VALENCE 1
  USEION K READ eK WRITE iK VALENCE 1
  USEION KCa WRITE iKCa VALENCE 1
  NONSPECIFIC_CURRENT iLeak
  RANGE gKbar, gK, gKCabar, gKCa
  RANGE gNabar, gNa
  RANGE gCaNbar, gCaN
  RANGE gLeakbar, eLeak
}

PARAMETER {
  alpha = 0.9 (/C/mm)
  kCa = 2 (/ms)
  gKCabar = 0.00003 (S/cm2)
  gNabar = 0.01 (S/cm2)
  gCaNbar = 0.0001 (S/cm2)
  gKbar = 0.0003 (S/cm2)
  gLeakbar = 0.00051 (S/cm2)
  eLeak = -65 (mV)
}

STATE {
  Ca  mNa  hNa  mCaN  hCaN  mK 
}

ASSIGNED {
  v (mV)
  eCa (mV)
  iCa (mA/cm2)
  eNa (mV)
  iNa (mA/cm2)
  eK (mV)
  iK (mA/cm2)
  iKCa (mA/cm2)
  gKCa (S/cm2)
  gK (S/cm2)
  gNa (S/cm2)
  gCaN (S/cm2)
  oldt (ms)
  oldv (mV)
  mNainf
  mNatau (ms)
  hNainf
  hNatau (ms)
  mCaNinf
  mCaNtau (ms)
  hCaNinf
  hCaNtau (ms)
  mKinf
  mKtau (ms)
  iLeak (mA/cm2)
}

BREAKPOINT {
  SOLVE states2 METHOD cnexp

  gNa = gNabar*mNa*mNa*mNa*hNa
  iNa = gNa*(v-eNa)
  gCaN = gCaNbar*mCaN*mCaN*hCaN
  iCa = gCaN*(v-eCa)
  gK = gKbar*mK*mK*mK*mK
  iK = gK*(v-eK)
  gKCa = gKCabar*Ca/(Ca+0.2)
  iKCa = gKCabar*(v-eK)
  iLeak = gLeakbar*(v-eLeak)
}

INITIAL {
  rates3(v)

  Ca = -alpha*iCa/kCa
  mNa = mNainf
  hNa = hNainf
  mCaN = mCaNinf
  hCaN = hCaNinf
  mK = mKinf
  eNa = 55
  eCa = 80
  eK = -80
}

DERIVATIVE states2 {
  rates3(v)

  Ca' = 0.01(-alpha*iCa-kCa*Ca)
  mNa' = (mNainf-mNa)/mNatau
  hNa' = (hNainf-hNa)/hNatau
  mCaN' = (mCaNinf-mCaN)/mCaNtau
  hCaN' = (hCaNinf-hCaN)/hCaNtau
  mK' = (mKinf-mK)/mKtau
}

UNITSOFF

PROCEDURE rates3(v(mV)) {
  TABLE mNainf, mNatau, hNainf, hNatau, mCaNinf, mCaNtau, 
   hCaNinf, hCaNtau, mKinf, mKtau FROM -100 TO 100 WITH 32000

  mNainf = NaPvtrapB(-(v+35),7.8)
  mNatau = 1e-2
  hNainf = NaPvtrapB((v+55),7)
  hNatau = 30/(NaPeFun((v+50),15)+NaPeFun(-(v+50),16))
  mCaNinf = NaPvtrapB(-(v+30),5)
  mCaNtau = 4
  hCaNinf = NaPvtrapB((v+45),5)
  hCaNtau = 40
  mKinf = NaPvtrapB(-(v+28),15)
  mKtau = 7/(NaPeFun((v+40),40)+NaPeFun(-(v+40),50))
}

FUNCTION NaPvtrapB(x,y) {
  LOCAL a
  a = x/y
  if (fabs(a)<1e-4) {
      NaPvtrapB = 1/(2+a)
  }else{
      NaPvtrapB = 1/(exp(a)+1)
  }
}

FUNCTION NaPeFun(x,y) {
  LOCAL a
  a=x/y
  if (fabs(a)<1e-4) {
      NaPeFun = 1+a
  }else{
      NaPeFun = exp(a)
  }
}

UNITSON